Evolutionary rescue experiment: effect of genetic diversity
Laure Olazcuaga
28-April-2022
1 DATA
1.1 Pop sive over time dataset
data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
head(data)## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5
## 1 no
## 2 no
## 3 no
## 4 no
## 5 no
## 6 notail(data)## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Med 5 6 8/25/21
## 572 Med 5 6 8/25/21
## 573 Med 5 6 8/25/21
## 574 Med 5 6 8/25/21
## 575 Med 5 6 8/25/21
## 576 Med 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5
## 571 no no yes
## 572 extinct yes yes
## 573 no no no
## 574 extinct yes yes
## 575 no no no
## 576 no no nodim(data)## [1] 576 23summary(data)## Block Old_Label Label Population
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## Length:576 Min. :0.0 Min. :1.0 Length:576
## Class :character 1st Qu.:1.0 1st Qu.:2.0 Class :character
## Mode :character Median :2.5 Median :3.5 Mode :character
## Mean :2.5 Mean :3.5
## 3rd Qu.:4.0 3rd Qu.:5.0
## Max. :5.0 Max. :6.0
##
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## MC_Box1 MC_Box2 MC_Total_Adults HC_Box1
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 26.52 1st Qu.: 25.99 1st Qu.: 51.05 1st Qu.: 22.00
## Median : 58.87 Median : 54.29 Median :101.50 Median : 56.50
## Mean : 75.88 Mean : 70.57 Mean :128.37 Mean : 72.06
## 3rd Qu.:106.93 3rd Qu.: 99.40 3rd Qu.:170.25 3rd Qu.:101.00
## Max. :327.40 Max. :299.00 Max. :514.40 Max. :288.00
## NA's :142 NA's :141 NA's :5 NA's :144
## HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## Min. : 0.00 Min. : 0.0 Min. : 0.00 Min. :0.0000
## 1st Qu.: 25.00 1st Qu.: 62.0 1st Qu.: 67.25 1st Qu.:0.4183
## Median : 52.00 Median :100.0 Median :100.00 Median :0.9386
## Mean : 68.08 Mean :132.1 Mean :138.30 Mean :1.4482
## 3rd Qu.: 96.75 3rd Qu.:169.0 3rd Qu.:188.00 3rd Qu.:2.3587
## Max. :290.00 Max. :499.0 Max. :499.00 Max. :6.8571
## NA's :142 NA's :44 NA's :122 NA's :142
## First_Throw Extinction_finalstatus Less_than_5
## Length:576 Length:576 Length:576
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## #Remove populations killed due to the pathogen
data_status <- data.frame(data$Population,data$Extinction_finalstatus)
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)## [1] 72 23data<-data[data$Extinction_finalstatus!="kill",]
#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove
#Check variables
str(data)## 'data.frame': 504 obs. of 26 variables:
## $ Block : chr "Block3" "Block3" "Block3" "Block3" ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : chr "High1" "High2" "High3" "High4" ...
## $ Genetic_Diversity : chr "High" "High" "High" "High" ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: chr "no" "no" "no" "no" ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)
#Order levels
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)## [1] "High" "Medium" "Low"data<-droplevels(data) #remove previous levels
str(data)## 'data.frame': 504 obs. of 26 variables:
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
## $ Genetic_Diversity : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...head(data)## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5 Nb_adults Lambda obs
## 1 no 100 NA 1
## 2 no 100 NA 2
## 3 no 100 NA 3
## 4 no 100 NA 4
## 5 no 100 NA 5
## 6 no 100 NA 6tail(data)## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Medium 5 6 8/25/21
## 572 Medium 5 6 8/25/21
## 573 Medium 5 6 8/25/21
## 574 Medium 5 6 8/25/21
## 575 Medium 5 6 8/25/21
## 576 Medium 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 571 no no yes 4 0.500000 499
## 572 extinct yes yes NA NA 500
## 573 no no no 46 1.314286 501
## 574 extinct yes yes NA NA 502
## 575 no no no 56 1.400000 503
## 576 no no no 133 3.243902 504dim(data)## [1] 504 261.2 Heterozygosity remain over time
#Upload He dataset
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]
#Merge dataset: add heterozygosity data to oridinal data
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)
### Creation of 4 He:
# 1- He_start: He was remaining for each population when we started the experiment
# 2- He_lost_gen_t: He was lost only during this specfic generation during the experiment (considering He=1 before this generation)
# 3- He_remain: He was remaining after each generation (He_remain[Gen=6]=He_end)
# 4- He_exp: He was lost during the entire experiment (considering He=1 at the beginning of the exp)
# 5- He_end: He was remaining for each population when we finished the experiment (considering the real He at the beginning of the exp, He_start, not 1)
#####
##### 1- He_start
#####
colnames(data)[27] <- "He_start"
#####
##### 2- He_lost_gen_t
#####
#Add He lost each generation
data$He_lost_gen_t <- 1 - (1/(2*data$Nb_adults))
#To consider extinct populations
is.na(data$He_lost_gen_t) <- sapply (data$He_lost_gen_t, is.infinite)
#####
##### 3- He_remain
#####
#Create new variable He remain per generation
data$He_remain <- NA
#Initiate with gen=1: He_start * He_lost_gen1
data$He_remain[data$Generation_Eggs=="1"] <- data$He_start[data$Generation_Eggs=="1"]*
data$He_lost_gen_t[data$Generation_Eggs=="1"]
#Compute for all the other genrations, except the first one
for(i in levels(data$Population)){
ifelse(sum(data$Generation_Eggs[data$Population==i])!=21,print("ERROR"),print(i))
for(j in 2:max(data$Generation_Eggs)){
data$He_remain[data$Population==i&data$Generation_Eggs==j] <-
data$He_remain[data$Population==i&data$Generation_Eggs==j-1]*
data$He_lost_gen_t[data$Population==i&data$Generation_Eggs==j]
}
}## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"#####
##### 4- He_exp
#####
# He at the end of the experiment
data$He_exp <- NA
#Compute the HE at the end of the experiment
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Eggs)!=21,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_adults))
#To consider extinct add this row:
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
#Compute for the whole experiment
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data$He_exp[data$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"#####
##### 5- He_end
#####
# Remaining He at the end of the experiment
data$He_end <- data$He_start * data$He_exp
#Save these values for Phenotyping
data_he <- merge(x = data_he,
y = data[data$Generation_Eggs==1,
c("Population","He_start",
"He_exp","He_end")],
by = "Population", all.x=TRUE)1.3 Phenotyping dataset
#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)
#Factors variables
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <- plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity,
levels = c("High", "Medium", "Low"))
#Add sum total
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx
data_phenotyping_all <- data_phenotyping
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]
#Merge with He dataset
data_phenotyping <- merge(x = data_phenotyping,
y = data_he, by = "Population", all.x=TRUE)
#Clean if less than 30 parents
pop_possible<-unique(data$Population[data$Generation_Eggs==5&data$Nb_adults>=30])
data_phenotyping <- data_phenotyping[data_phenotyping$Population %in% pop_possible,]1.4 Old He remain
#Upload growth rate end of experiment
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]
### Additional: He remaining
# New vector for the lost during exp
data_he$He_remain_exp <- NA
# He lost during exp
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))
#To consider extinct add this row:
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"# Remaining He at the end of the experiment
data_he$He_end <- data_he$He_remain * data_he$He_remain_exp
# Remove kill population
data_he <- data_he[!is.na(data_he$He_remain_exp),]
data_he <- droplevels(data_he)
## Merge dataset He
data_old_merged <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)2 PLOT
2.1 Population size
PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Population,
colour = Genetic_Diversity)) +
geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
SUM_Popsize <- Rmisc::summarySE(data,
measurevar = c("Nb_adults"),
groupvars = c("Generation_Eggs",
"Genetic_Diversity"),
na.rm=TRUE)
SUM_Popsize## Generation_Eggs Genetic_Diversity N Nb_adults sd se ci
## 1 1 High 27 100.00000 0.00000 0.000000 0.00000
## 2 1 Medium 23 100.00000 0.00000 0.000000 0.00000
## 3 1 Low 34 100.00000 0.00000 0.000000 0.00000
## 4 2 High 27 344.29630 73.77294 14.197610 29.18360
## 5 2 Medium 23 282.04348 100.75735 21.009360 43.57075
## 6 2 Low 34 282.61765 78.53006 13.467794 27.40043
## 7 3 High 27 151.85185 80.96899 15.582489 32.03027
## 8 3 Medium 23 118.73913 88.54491 18.462891 38.28969
## 9 3 Low 34 104.61765 85.13341 14.600260 29.70445
## 10 4 High 26 114.00000 80.20224 15.728954 32.39439
## 11 4 Medium 23 62.65217 64.50483 13.450188 27.89398
## 12 4 Low 34 46.38235 39.63241 6.796903 13.82840
## 13 5 High 27 103.62963 51.56486 9.923661 20.39838
## 14 5 Medium 20 61.50000 50.59904 11.314290 23.68108
## 15 5 Low 31 51.51613 46.13377 8.285870 16.92200
## 16 6 High 27 147.66667 42.60282 8.198916 16.85311
## 17 6 Medium 19 69.73684 46.02396 10.558620 22.18284
## 18 6 Low 28 77.03571 46.14428 8.720450 17.89289# SUM_Popsize <- Rmisc::summarySE(data[data$Extinction_finalstatus=="no",],
# measurevar = c("Nb_adults"),
# groupvars = c("Generation_Eggs",
# "Genetic_Diversity"),
# na.rm=TRUE)
# SUM_Popsize
PLOT_Pop_size_average <- PLOT_Pop_size +
geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 5, position = position_dodge(0.4)) +
geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
ymin = Nb_adults-ci, ymax = Nb_adults+ci,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 1, width=.2, position = position_dodge(0.4)) +
geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
linetype = "dashed", size = 1.5, position = position_dodge(0.4))
PLOT_Pop_size_average
# cowplot::save_plot(file = here::here("figures","PopulationSize_average.pdf"), PLOT_Pop_size_average, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)
#Prediction with models
# Test effect
data$gen_square <- data$Generation_Eggs^2
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2&data$Extinction_finalstatus=="no",],
family = "poisson")
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(2,6, length = 100),each=3),
Block = "Block4",
gen_square = (rep(seq(2,6, length = 100),each=3))^2,
Estimates = NA)
filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)
#Predcit estimate minimun values
tapply(filldata$Estimates, filldata$Genetic_Diversity, min)## High Low Medium
## 85.85198 44.16308 51.99642filldata[filldata$Estimates>=44.16&filldata$Estimates<=44.17,]## Genetic_Diversity Generation_Eggs Block gen_square Estimates
## 198 Low 4.626263 Block4 21.40231 44.16308filldata[filldata$Estimates>=51.99&filldata$Estimates<=52.00,]## Genetic_Diversity Generation_Eggs Block gen_square Estimates
## 230 Medium 5.070707 Block4 25.71207 51.99642filldata[filldata$Estimates>=85.85&filldata$Estimates<=85.86,]## Genetic_Diversity Generation_Eggs Block gen_square Estimates
## 181 High 4.424242 Block4 19.57392 85.85198## Change name vector
vector_names <- c(`Low` = "Extreme-bottleneck",
`Medium` = "Severe-bottleneck",
`High` = "Diverse-source")
PLOT_Pop_size_predict<-ggplot2::ggplot(data[data$Extinction_finalstatus=="no",],
aes(x = factor(Generation_Eggs),
y = Nb_adults)) +
geom_point(aes(group = Population,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 0.7, alpha = 0.5) +
geom_line(aes(group = Population,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 0.05, alpha = 0.5) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.4) +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse-source","Severe-bottleneck","Extreme-bottleneck"),
values = c("#00AFBB", "#E7B800", "#FC4E07")) +
ylab("Population size") +
xlab("Generation") +
theme_LO
PLOT_Pop_size_predict
#
# cowplot::save_plot(file = here::here("figures","PopulationSize_predict.pdf"), PLOT_Pop_size_predict, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)2.2 Extinction
#Percent of extinction
length(unique(data$Population[data$Genetic_Diversity=="Low"&
data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="Low"]))## [1] 0.2352941length(unique(data$Population[data$Genetic_Diversity=="Medium"&
data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="Medium"]))## [1] 0.2173913length(unique(data$Population[data$Genetic_Diversity=="High"&
data$Extinction_finalstatus=="yes"]))/
length(unique(data$Population[data$Genetic_Diversity=="High"]))## [1] 0#Create data with percent
vector_names <- c(`Low` = "Extreme-bottleneck",
`Medium` = "Severe-bottleneck",
`High` = "Diverse-source")
f_labels = data.frame(Genetic_Diversity = c("High","Medium","Low"),
label = c("Extinction = 0 %", "Extinction = 21.7 %", "Extinction = 23.5 %"))
f_labels$Genetic_Diversity <- factor(f_labels$Genetic_Diversity,
levels = c("High", "Medium", "Low"))
## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults)) +
facet_wrap(~ Genetic_Diversity, ncol=3, labeller = ggplot2::as_labeller(vector_names)) +
#geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
geom_line(aes(group = Population,
colour = Extinction_finalstatus),
position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
geom_text(x = 5, y = 420,
aes(label = label),
data = f_labels,
col="red",
size = 3,
vjust = 0) +
scale_color_manual(values = c("black", "red")) +
ylab("Population size") +
xlab("Generation") +
theme(legend.position = "none") +
theme_LO
PLOT_Extinction
cowplot::save_plot(here::here("figures","Extinction.pdf"), PLOT_Extinction, base_height = 8/cm(1),
base_width = 20/cm(1), dpi = 610)
#
PLOT_Extinction
2.3 Growth rate G2
tapply(data[data$Generation_Eggs==2,]$Lambda,data[data$Generation_Eggs==2,]$Genetic_Diversity,mean)## High Medium Low
## 3.442963 2.820435 2.826176PLOT_Growth<-ggplot2::ggplot(data[data$Generation_Eggs==2,], aes(x = Genetic_Diversity, y = Lambda,
colour = Genetic_Diversity)) +
geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Genetic diversity") +
theme_LO
PLOT_Growth
#
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"), PLOT_Growth, base_height = 10/cm(1),
# base_width = 8/cm(1), dpi = 610)2.4 Life stage proportion
SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_adults"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_adults## Genetic_Diversity Week N p_adults sd se ci
## 1 High 4-week 50 0.6508098 0.17776618 0.025139934 0.05052059
## 2 High 5-week 105 0.9632040 0.05259208 0.005132461 0.01017786
## 3 Medium 4-week 24 0.7121270 0.15722597 0.032093616 0.06639070
## 4 Medium 5-week 49 0.9563136 0.06160975 0.008801393 0.01769639
## 5 Low 4-week 30 0.6568606 0.19029917 0.034743716 0.07105888
## 6 Low 5-week 59 0.9582021 0.07845154 0.010213520 0.02044458SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_pupae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_pupae## Genetic_Diversity Week N p_pupae sd se ci
## 1 High 4-week 50 0.28825555 0.15629878 0.022103985 0.044419621
## 2 High 5-week 105 0.01458769 0.02407445 0.002349426 0.004658999
## 3 Medium 4-week 24 0.23019148 0.14506173 0.029610601 0.061254195
## 4 Medium 5-week 49 0.02283167 0.03485130 0.004978757 0.010010462
## 5 Low 4-week 30 0.29380227 0.17562943 0.032065401 0.065581108
## 6 Low 5-week 59 0.01675308 0.02459640 0.003202178 0.006409856SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_larvae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_larvae## Genetic_Diversity Week N p_larvae sd se ci
## 1 High 4-week 50 0.06093466 0.04274850 0.006045551 0.012148989
## 2 High 5-week 105 0.02220828 0.04324043 0.004219833 0.008368088
## 3 Medium 4-week 24 0.05768150 0.05999020 0.012245449 0.025331640
## 4 Medium 5-week 49 0.02085474 0.04462923 0.006375605 0.012819012
## 5 Low 4-week 30 0.04933713 0.08198954 0.014969174 0.030615398
## 6 Low 5-week 59 0.02504479 0.06691319 0.008711355 0.017437671SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
p_pupae=SUM_Prop_pupae$p_pupae,
p_larvae=SUM_Prop_larvae$p_larvae)
SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop## Genetic_Diversity Week N Stage Proportion
## 1 High 4-week 50 p_adults 0.65080979
## 2 High 5-week 105 p_adults 0.96320403
## 3 Medium 4-week 24 p_adults 0.71212702
## 4 Medium 5-week 49 p_adults 0.95631359
## 5 Low 4-week 30 p_adults 0.65686059
## 6 Low 5-week 59 p_adults 0.95820213
## 7 High 4-week 50 p_pupae 0.28825555
## 8 High 5-week 105 p_pupae 0.01458769
## 9 Medium 4-week 24 p_pupae 0.23019148
## 10 Medium 5-week 49 p_pupae 0.02283167
## 11 Low 4-week 30 p_pupae 0.29380227
## 12 Low 5-week 59 p_pupae 0.01675308
## 13 High 4-week 50 p_larvae 0.06093466
## 14 High 5-week 105 p_larvae 0.02220828
## 15 Medium 4-week 24 p_larvae 0.05768150
## 16 Medium 5-week 49 p_larvae 0.02085474
## 17 Low 4-week 30 p_larvae 0.04933713
## 18 Low 5-week 59 p_larvae 0.02504479SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))
PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity,
y = Proportion, fill = Stage)) +
facet_wrap(~ Week, ncol=2) +
geom_bar(stat="identity") +
xlab("Genetic history") +
labs(color="Stage of individuals") +
ylab("Proportion of each stage") +
scale_fill_manual(values= c("#619CFF","#F8766D","#00BA38"),
breaks = c("p_adults", "p_pupae", "p_larvae"),
labels = c( "Adult", "Pupae","Larvae")) +
ggplot2::scale_x_discrete(breaks=c("High", "Medium", "Low"),
labels=c("Diverse-\nsource",
"Severe-\nbottleneck",
"Extreme-\nbottleneck")) +
theme_LO
PLOT_prop
# cowplot::save_plot(here::here("figures","Life_stage_proportion.pdf"), PLOT_prop,
# base_height = 10/cm(1), base_width = 18/cm(1), dpi = 610)2.5 Remaining heterozygosity
# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_phenotyping,
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity",
"He_end",
"Population"),
na.rm=TRUE)
PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_end,
y = Lambda,
colour = Genetic_Diversity,
size = N)) +
geom_point(alpha = 0.8) +
ylab("Growth rate") +
xlab("Expected heterozygosity") +
scale_color_manual(name="Genetic history",
breaks = c("High", "Medium", "Low"),
labels = c("Diverse-source","Severe-bottleneck","Extreme-bottleneck"),
values = c("#00AFBB", "#E7B800", "#FC4E07")) +
labs(size="Replicates") +
theme_LO
PLOT_He_Final
#
#
# cowplot::save_plot(here::here("figures","Heterozygosity_Growthrate.pdf"),
# PLOT_He_Final,
# base_height = 12/cm(1), base_width = 16/cm(1), dpi = 610)
#
# # ALL WITHOUT PSEUDO
# SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
# measurevar = c("Lambda"),
# groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_end",
# "Population"),
# na.rm=TRUE)
#
# SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_end)
#
# SUM_POP_He_ALL <- SUM_POP_He_ALL[!is.na(SUM_POP_He_ALL$He_end),]
# SUM_POP_He_ALL <- SUM_POP_He_ALL[SUM_POP_He_ALL$HE!="-Inf",]
#
#
# facet_names <- c('Initial'="Before rescue experiment",'Final'="After rescue experiment")
#
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
# colour = Genetic_Diversity,
# shape = Phenotyping)) +
# #facet_wrap(~ Phenotyping, ncol=1) +
# geom_point(size = 3, alpha = 0.8) +
# scale_shape_manual(values = c(1, 19)) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL
#
#
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
# colour = Genetic_Diversity,
# shape = Phenotyping)) +
# facet_wrap(~ Phenotyping, ncol=1) +
# geom_point(size = 3, alpha = 0.8) +
# scale_shape_manual(values = c(1, 19)) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL
#
#
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
# colour = Genetic_Diversity)) +
# facet_wrap(~ Phenotyping, ncol=1) +
# geom_point(size = 3, alpha = 0.8) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL2.6 He over time
#Compute for all the other genrations, except the first one
data$ID_extinction <- "extant"
for(i in unique(data$Population[data$Extinction_finalstatus=="yes"])){
gen_all <- data$Generation_Eggs[data$Population==i&!is.na(data$He_remain)]
maxgen <- max(gen_all)
data$ID_extinction[data$Population==i&data$Generation_Eggs==maxgen] <- "willextinct"
}
data[data$ID_extinction=="willextinct",]## Population Block Old_Label Label
## 205 Low16 Block4 B4-O1 7/14 BE B4 | G5 Low 16
## 273 Low27 Block5 B5-B1 6/16 BE B5 | G4 Low 27
## 278 Low28 Block5 B5_E1 5/12 BE B5 | G3 Low 28
## 299 Low30 Block5 B5-D1 6/16 BE B5 | G4 Low 30
## 303 Low31 Block5 B(2)5-P1 6/16 BE B5 | G4 Low 31
## 310 Low32 Block5 B(2)5_Q1 5/12 BE B5 | G3 Low 32
## 317 Low33 Block5 B(2)5-T1 7/21 BE B5 | G5 Low 33
## 336 Low36 Block5 B(2)5-X1 6/16 BE B5 | G4 Low 36
## 402 Med14 Block4 B4-Backup Fam L 6/9 BE B4 | G4 Med 14
## 409 Med17 Block5 B5_Backup_Fam_F 5/12 BE B5 | G3 Med 17
## 437 Med20 Block5 B5 - Backup Fam I 6/16 BE B5 | G4 Med 20
## 449 Med22 Block5 B5_Backup_Fam_B 5/12 BE B5 | G3 Med 22
## 479 Med5 Block3 B3 - Backup Fam E 7/7 BE B3 | G5 Med 5
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 205 Low 4 5 7/14/21
## 273 Low 3 4 6/16/21
## 278 Low 2 3 5/12/21
## 299 Low 3 4 6/16/21
## 303 Low 3 4 6/16/21
## 310 Low 2 3 5/12/21
## 317 Low 4 5 7/21/21
## 336 Low 3 4 6/16/21
## 402 Medium 3 4 6/9/21
## 409 Medium 2 3 5/12/21
## 437 Medium 3 4 6/16/21
## 449 Medium 2 3 5/12/21
## 479 Medium 4 5 7/7/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 205 7/14/21 7/15/21 <NA> DSC_0994 0.000000 2.000000
## 273 6/16/21 6/17/21 DSC_0526 DSC_0527 5.090878 0.000000
## 278 5/12/21 5/13/21 DSC_0918 DSC_0919 3.085768 3.638847
## 299 6/16/21 6/17/21 DSC_0559 <NA> 1.383920 0.000000
## 303 6/16/21 6/17/21 <NA> DSC_0541 0.000000 1.000000
## 310 5/12/21 5/13/21 DSC_0953 DSC_0954 88.254170 70.952124
## 317 7/21/21 7/22/21 DSC_0117 DSC_0118 1.000000 1.000000
## 336 6/16/21 6/17/21 DSC_0540 <NA> 1.000000 0.000000
## 402 6/9/21 6/10/21 DSC_0376 DSC_0377 10.974990 1.000000
## 409 5/12/21 5/13/21 DSC_0916 DSC_0917 1.543122 1.084327
## 437 6/16/21 6/17/21 DSC_0542 <NA> 1.000000 0.000000
## 449 5/12/21 5/13/21 DSC_0922 DSC_0923 9.168447 6.251069
## 479 7/7/21 7/8/21 <NA> DSC_0830 NA NA
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 205 2.0 0 2 2 5 0.400000000
## 273 5.1 3 0 3 19 0.157894737
## 278 6.7 1 1 2 374 0.005347594
## 299 1.4 1 0 1 146 0.006849315
## 303 1.0 1 1 2 272 0.007352941
## 310 159.2 71 56 127 276 0.460144928
## 317 2.0 1 1 2 3 0.666666667
## 336 1.0 1 0 1 22 0.045454545
## 402 12.0 7 1 8 138 0.057971014
## 409 2.6 3 2 5 416 0.012019231
## 437 1.0 1 0 1 20 0.050000000
## 449 15.4 10 6 16 200 0.080000000
## 479 0.0 NA NA 1 11 0.090909091
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 205 no yes yes 2 0.400000000 378
## 273 no yes yes 3 0.157894737 319
## 278 no yes yes 2 0.005347594 236
## 299 no yes yes 1 0.006849315 322
## 303 no yes yes 2 0.007352941 323
## 310 no yes yes 127 0.460144928 240
## 317 no yes yes 2 0.666666667 409
## 336 no yes yes 1 0.045454545 328
## 402 no yes yes 8 0.057971014 308
## 409 no yes yes 5 0.012019231 245
## 437 no yes yes 1 0.050000000 332
## 449 no yes yes 16 0.080000000 250
## 479 no yes yes 1 0.090909091 359
## He_start He_lost_gen_t He_remain He_exp He_end gen_square
## 205 0.5544299 0.7500000 0.3575672 0.6449276 0.3575672 25
## 273 0.5367998 0.8333333 0.4324691 0.8056432 0.4324691 16
## 278 0.5386585 0.7500000 0.4014365 0.7452523 0.4014365 9
## 299 0.5540444 0.5000000 0.2742819 0.4950540 0.2742819 16
## 303 0.5532775 0.7500000 0.4116634 0.7440450 0.4116634 16
## 310 0.5528880 0.9960630 0.5469651 0.9892872 0.5469651 9
## 317 0.5523333 0.7500000 0.3359893 0.6083089 0.3359893 25
## 336 0.5427371 0.5000000 0.2635269 0.4855518 0.2635269 16
## 402 0.6802331 0.9375000 0.6315974 0.9285014 0.6315974 16
## 409 0.6825509 0.9000000 0.6104897 0.8944237 0.6104897 9
## 437 0.6819650 0.5000000 0.3302071 0.4841994 0.3302071 16
## 449 0.6809168 0.9687500 0.6546991 0.9614965 0.6546991 9
## 479 0.6807065 0.5000000 0.3209456 0.4714889 0.3209456 25
## ID_extinction
## 205 willextinct
## 273 willextinct
## 278 willextinct
## 299 willextinct
## 303 willextinct
## 310 willextinct
## 317 willextinct
## 336 willextinct
## 402 willextinct
## 409 willextinct
## 437 willextinct
## 449 willextinct
## 479 willextinctdim(data[data$ID_extinction=="willextinct",])## [1] 13 33data[data$ID_extinction=="willextinct","Population"]## [1] Low16 Low27 Low28 Low30 Low31 Low32 Low33 Low36 Med14 Med17 Med20 Med22
## [13] Med5
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9data$Population[data$Genetic_Diversity=="Medium"&data$Extinction_finalstatus=="yes"]## [1] Med14 Med14 Med14 Med14 Med14 Med14 Med17 Med17 Med17 Med17 Med17 Med17
## [13] Med20 Med20 Med20 Med20 Med20 Med20 Med22 Med22 Med22 Med22 Med22 Med22
## [25] Med5 Med5 Med5 Med5 Med5 Med5
## 84 Levels: High1 High13 High14 High15 High16 High17 High18 High19 ... Med9PLOT_He_extinction <-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = He_remain,
group = Population,
colour = Extinction_finalstatus,
shape = ID_extinction)) +
geom_point(position = position_dodge(0.5), size = 3, alpha = 0.6) +
geom_line(position = position_dodge(0.5), size = 0.25, alpha = 0.4) +
ylab("Expected heterozygozity") +
scale_color_manual(values = c("black", "red")) +
scale_shape_manual(values = c(16, 13), guide=FALSE) +
labs(color="Extinct populations") +
xlab("Generation") +
theme_LO
PLOT_He_extinction
#
# cowplot::save_plot(here::here("figures","Heterozygosity_over_time.pdf"),
# PLOT_He_extinction,
# base_height = 12/cm(1), base_width = 14/cm(1), dpi = 610)3 ANALYSIS
3.1 Probability of extinction
############
###### Clean dataset
############
#Prepare clean dataset
data_proba_extinction <- data[data$Generation_Eggs==6,c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)
############
###### Analysis
############
#Analysis
mod0 <- glm(y ~ Genetic_Diversity + Block,
data = data_proba_extinction, family = binomial)
summary(mod0)##
## Call:
## glm(formula = y ~ Genetic_Diversity + Block, family = binomial,
## data = data_proba_extinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26011 -0.44258 -0.30957 -0.00005 2.47474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.3144 1855.7435 -0.011 0.99084
## Genetic_DiversityMedium 18.3002 1855.7433 0.010 0.99213
## Genetic_DiversityLow 18.5062 1855.7432 0.010 0.99204
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 72.388 on 83 degrees of freedom
## Residual deviance: 46.833 on 79 degrees of freedom
## AIC: 56.833
##
## Number of Fisher Scoring iterations: 18#Test genetic diversity effect
mod1 <- glm(y ~ Block ,
data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 81 58.903
## 2 79 46.833 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#We keep genetic diversity
# Perform the lrtest
(A <- logLik(mod1))## 'log Lik.' -29.45146 (df=3)(B <- logLik(mod0))## 'log Lik.' -23.41669 (df=5)#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))## [1] 12.06954#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))## [1] 0.002394043lmtest::lrtest(mod1, mod0)## Likelihood ratio test
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -29.451
## 2 5 -23.417 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,
data = data_proba_extinction, family = binomial)
summary(mod0)##
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial,
## data = data_proba_extinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26011 -0.44258 -0.30957 -0.00005 2.47474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Genetic_DiversityHigh -21.3144 1855.7435 -0.011 0.99084
## Genetic_DiversityMedium -3.0142 1.1293 -2.669 0.00760 **
## Genetic_DiversityLow -2.8082 1.0659 -2.635 0.00843 **
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.449 on 84 degrees of freedom
## Residual deviance: 46.833 on 79 degrees of freedom
## AIC: 56.833
##
## Number of Fisher Scoring iterations: 18############
###### Predicted value
############
#Extract probability of extinction
(p_high <- plogis(-21.3144))## [1] 5.536989e-10(p_med <- plogis(-3.0142))## [1] 0.04678847(p_low <- plogis(-2.8082))## [1] 0.05688267data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
Block = "Block4")
data_predict_extinct$predict <- predict(mod0, type = "response",
re.form = NA,
newdata = data_predict_extinct)
data_predict_extinct## Genetic_Diversity Block predict
## 1 High Block4 1.160723e-09
## 2 Medium Block4 9.329490e-02
## 3 Low Block4 1.122446e-01p_high## [1] 5.536989e-10p_med## [1] 0.04678847p_low## [1] 0.05688267############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High -20.07 1855.743 Inf -4451.10 4410.961
## Medium -1.77 0.650 Inf -3.32 -0.216
## Low -1.56 0.532 Inf -2.83 -0.290
##
## Results are averaged over the levels of: Block
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium -18.300 1855.743 Inf -0.010 0.9999
## High - Low -18.506 1855.743 Inf -0.010 0.9999
## Medium - Low -0.206 0.751 Inf -0.274 0.9594
##
## Results are averaged over the levels of: Block
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.2 Population size G6
############
###### Analysis
############
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),
# data = data[data$Generation_Eggs==6,], family = "poisson")
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block,
data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])
mod1 <- lm(log(Nb_adults) ~ Block ,
data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])
anova(mod0, mod1) # Effect of genetic diversity ## Analysis of Variance Table
##
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 66 29.706
## 2 68 42.118 -2 -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block,
data = data[data$Generation_Eggs==6&data$Extinction_finalstatus == "no",])
summary(mod0)##
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data[data$Generation_Eggs ==
## 6 & data$Extinction_finalstatus == "no", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2417 -0.2947 0.1088 0.4381 1.2624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.11532 0.17980 28.451 < 2e-16 ***
## Genetic_DiversityMedium -0.94813 0.20540 -4.616 1.86e-05 ***
## Genetic_DiversityLow -0.80532 0.18719 -4.302 5.72e-05 ***
## BlockBlock4 -0.02077 0.18447 -0.113 0.9107
## BlockBlock5 -0.53922 0.21414 -2.518 0.0142 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared: 0.3362, Adjusted R-squared: 0.2959
## F-statistic: 8.356 on 4 and 66 DF, p-value: 1.623e-05############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 4.93 0.130 66 4.61 5.25
## Medium 3.98 0.159 66 3.59 4.37
## Low 4.12 0.136 66 3.79 4.46
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.948 0.205 66 4.616 0.0001
## High - Low 0.805 0.187 66 4.302 0.0002
## Medium - Low -0.143 0.207 66 -0.690 0.7704
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.3 Population size over time
#######################
####################### MODEL INCLUDING GENETIC DIVERSITY WITHOUT GEN 1
#######################
###################### REMOVE GEN 1 ##########################33
#If we dont consider Gen=1
PLOT_Pop_size_average
fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2,])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
fit <- glm(Nb_adults ~ Generation_Eggs,
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit6 <- glm(Nb_adults~exp(Generation_Eggs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit7 <- glm(Nb_adults~log(Generation_Eggs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
#generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(1,6, length=100)
# plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
# log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
# lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
# lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
# lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
# lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
# lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
# lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
data[data$Generation_Eggs>=2,]$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")
# Check the r squared
rsq::rsq(fit,adj=TRUE)## [1] 0.4442203rsq::rsq(fit2,adj=TRUE)## [1] 0.5938535rsq::rsq(fit3,adj=TRUE)## [1] 0.5928435rsq::rsq(fit4,adj=TRUE)## [1] 0.592052rsq::rsq(fit5,adj=TRUE)## [1] 0.592052rsq::rsq(fit6,adj=TRUE)## [1] 0.121301rsq::rsq(fit7,adj=TRUE)## [1] 0.5205104# Best model is fit 2
data$gen_square <- data$Generation_Eggs^2
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~Generation_Eggs + gen_square,
data = data[data$Generation_Eggs>=2,], family = "poisson")
# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
#Add the interaction with Genetic diversity
fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
gen_square*Genetic_Diversity,
data = data, family = "poisson")
summary(fit2)##
## Call:
## glm(formula = Nb_adults ~ Generation_Eggs * Genetic_Diversity +
## gen_square * Genetic_Diversity, family = "poisson", data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -15.961 -6.480 -2.662 5.254 21.115
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.070332 0.026463 191.604 < 2e-16
## Generation_Eggs 0.156559 0.017811 8.790 < 2e-16
## Genetic_DiversityMedium -0.164697 0.042046 -3.917 8.96e-05
## Genetic_DiversityLow 0.024718 0.037632 0.657 0.511290
## gen_square -0.036751 0.002566 -14.324 < 2e-16
## Generation_Eggs:Genetic_DiversityMedium 0.080423 0.029501 2.726 0.006409
## Generation_Eggs:Genetic_DiversityLow -0.096063 0.026425 -3.635 0.000278
## Genetic_DiversityMedium:gen_square -0.035166 0.004448 -7.905 2.67e-15
## Genetic_DiversityLow:gen_square -0.009935 0.003956 -2.511 0.012038
##
## (Intercept) ***
## Generation_Eggs ***
## Genetic_DiversityMedium ***
## Genetic_DiversityLow
## gen_square ***
## Generation_Eggs:Genetic_DiversityMedium **
## Generation_Eggs:Genetic_DiversityLow ***
## Genetic_DiversityMedium:gen_square ***
## Genetic_DiversityLow:gen_square *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 39168 on 486 degrees of freedom
## Residual deviance: 30719 on 478 degrees of freedom
## (17 observations deleted due to missingness)
## AIC: 33747
##
## Number of Fisher Scoring iterations: 5####### Add genetic diversity
#Test oversdispersion
mod0 <- glm(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block ,
data = data[data$Generation_Eggs>=2,], family = "poisson")
sqrt(deviance(mod0)/(nobs(mod0)-length(coef(mod0))))## [1] 6.395475sigma(mod0)## [1] 6.395475# ccl: overdispersion
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
blmeco::dispersion_glmer(mod1)## [1] 1.076696sigma(mod1)## [1] 1#Convergence issue
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001## [1] 0.004551091# Test effect
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
anova(mod1, mod2, test ="Chisq") ## Data: data[data$Generation_Eggs >= 2, ]
## Models:
## mod2: Nb_adults ~ Generation_Eggs + gen_square + Genetic_Diversity +
## mod2: Block + (1 | obs)
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## mod1: Genetic_Diversity + Block + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod2 8 4669.1 4701.1 -2326.6 4653.1
## mod1 12 4660.8 4708.8 -2318.4 4636.8 16.376 4 0.002554 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(2,6, length = 100),each=3),
Block = "Block4",
gen_square = (rep(seq(2,6, length = 100),each=3))^2,
Estimates = NA)
filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)
#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data[data$Generation_Eggs>=2,]) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.3) +
geom_point(data = data,
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
emmeans::emmeans(mod1, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.88 0.0790 Inf 4.69 5.07
## Medium 4.19 0.0885 Inf 3.98 4.40
## Low 4.04 0.0733 Inf 3.87 4.22
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.688 0.119 Inf 5.786 <.0001
## High - Low 0.835 0.107 Inf 7.778 <.0001
## Medium - Low 0.147 0.115 Inf 1.287 0.4028
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates###################################################################
###################################################################
###################################################################
#################### INCLUDING GEN 1 ###########################
PLOT_Pop_size_average
fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")
#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")
# Check the r squared
summary(fit)$adj.r.squared## NULLsummary(fit2)$adj.r.squared## NULLsummary(fit3)$adj.r.squared## NULLsummary(fit4)$adj.r.squared # best r square## NULLsummary(fit5)$adj.r.squared ## NULLsummary(fit6)$adj.r.squared ## NULLsummary(fit7)$adj.r.squared ## NULLrsq::rsq(fit,adj=TRUE)## [1] 0.1067643rsq::rsq(fit2,adj=TRUE)## [1] 0.1486733rsq::rsq(fit3,adj=TRUE)## [1] 0.5240033rsq::rsq(fit4,adj=TRUE)## [1] 0.5980368rsq::rsq(fit5,adj=TRUE)## [1] 0.5989426rsq::rsq(fit6,adj=TRUE)## [1] 0.07185974rsq::rsq(fit7,adj=TRUE)## [1] 0.04937507#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#######################
####################### MODEL INCLUDING GENETIC DIVERSITY and GEN 1
#######################
#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)## [1] 5.760264mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)## [1] 1.083377#Convergence issue
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001## [1] 0.06384154#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(1,6, length = 100),each=3),
Block = "Block4",
Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")
#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.3) +
geom_point(data = data,
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity +
# gen_square*Genetic_Diversity,
# data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity +
# gen_square,
# data = data[data$Generation_Eggs>=2,])
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## mod1: (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## mod: (1 | Block) + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 9 5611.2 5648.9 -2796.6 5593.2
## mod 17 5594.8 5666.0 -2780.4 5560.8 32.47 8 7.672e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Signifcant interaction
summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 5594.8 5666.0 -2780.4 5560.8 470
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.03685 -0.03969 0.00800 0.05695 0.26801
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.64646 0.8040
## Block (Intercept) 0.07945 0.2819
## Number of obs: 487, groups: obs, 487; Block, 3
##
## Fixed effects:
## Estimate
## (Intercept) -1.993477
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 10.858800
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.147603
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.942314
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.058998
## Genetic_DiversityMedium -0.352704
## Genetic_DiversityLow -0.787854
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.772686
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.503761
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.095237
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.006018
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 1.489498
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow -0.833931
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.138807
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow -0.007028
## Std. Error
## (Intercept) 1.238244
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 1.985855
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 1.032426
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.213266
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.015171
## Genetic_DiversityMedium 1.833248
## Genetic_DiversityLow 1.630486
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 2.976126
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.553127
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.322073
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.022992
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 2.637956
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 1.373155
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.284239
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.020268
## z value
## (Intercept) -1.610
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 5.468
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -4.986
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 4.418
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -3.889
## Genetic_DiversityMedium -0.192
## Genetic_DiversityLow -0.483
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.260
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.324
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.296
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.262
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.565
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow -0.607
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.488
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow -0.347
## Pr(>|z|)
## (Intercept) 0.107415
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 4.55e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 6.17e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 9.94e-06
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.000101
## Genetic_DiversityMedium 0.847434
## Genetic_DiversityLow 0.628952
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.795151
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.745671
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.767459
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.793538
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.572318
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.543645
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.625306
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.728786
##
## (Intercept)
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.398806 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?summary(mod1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 5611.2 5648.9 -2796.6 5593.2 478
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05397 -0.05185 0.01426 0.07054 0.21832
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.69318 0.8326
## Block (Intercept) 0.07665 0.2769
## Number of obs: 487, groups: obs, 487; Block, 3
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.01996 0.81642 -2.474
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.74179 1.30992 8.964
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.66290 0.68789 -8.232
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 1.03332 0.14306 7.223
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.06412 0.01022 -6.274
## Genetic_DiversityMedium -0.55984 0.10012 -5.592
## Genetic_DiversityLow -0.67366 0.09048 -7.445
## Pr(>|z|)
## (Intercept) 0.0134 *
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.08e-13 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 3.52e-10 ***
## Genetic_DiversityMedium 2.25e-08 ***
## Genetic_DiversityLow 9.69e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.965
## s::(G_E,4,r=TRUE)2 0.941 -0.993
## s::(G_E,4,r=TRUE)3 -0.915 0.976 -0.995
## s::(G_E,4,r=TRUE)4 0.889 -0.956 0.984
## Gntc_DvrstM -0.062 0.007 -0.008
## Gntc_DvrstL -0.065 0.004 -0.005
## s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1
## s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)3
## s::(G_E,4,r=TRUE)4 -0.997
## Gntc_DvrstM 0.008 -0.009
## Gntc_DvrstL 0.005 -0.006 0.494
## convergence code: 0
## Model failed to converge with max|grad| = 0.053402 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.55 0.206 Inf 4.06 5.04
## Medium 3.93 0.216 Inf 3.42 4.45
## Low 3.69 0.201 Inf 3.21 4.17
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.619 0.187 Inf 3.310 0.0027
## High - Low 0.861 0.168 Inf 5.115 <.0001
## Medium - Low 0.242 0.184 Inf 1.319 0.3844
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#######################
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
#######################
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")
PLOT_Pop_size <- ggplot2::ggplot(data = data) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates_Withoutextpop,
colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
geom_point(data = data[data$Extinction_finalstatus=="no",],
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs),
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity +
(1|Block) + (1|obs),
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## mod1: (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## mod: (1 | Block) + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 9 4759.9 4796.3 -2370.9 4741.9
## mod 17 4756.3 4825.1 -2361.1 4722.3 19.599 8 0.01196 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data[data$Extinction_finalstatus == "no", ]
##
## AIC BIC logLik deviance df.resid
## 4756.3 4825.1 -2361.1 4722.3 407
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.57063 -0.04879 0.00393 0.06967 0.29650
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.36624 0.6052
## Block (Intercept) 0.02456 0.1567
## Number of obs: 424, groups: obs, 424; Block, 3
##
## Fixed effects:
## Estimate
## (Intercept) -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.0583188
## Genetic_DiversityMedium 1.1892488
## Genetic_DiversityLow 0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.0006044
## Std. Error
## (Intercept) 0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.0120893
## Genetic_DiversityMedium 1.5049889
## Genetic_DiversityLow 1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.0173582
## z value
## (Intercept) -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -4.824
## Genetic_DiversityMedium 0.790
## Genetic_DiversityLow 0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.035
## Pr(>|z|)
## (Intercept) 0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 1.41e-06
## Genetic_DiversityMedium 0.4294
## Genetic_DiversityLow 0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.9722
##
## (Intercept) *
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.53 0.133 Inf 4.21 4.85
## Medium 4.30 0.151 Inf 3.94 4.66
## Low 4.01 0.137 Inf 3.68 4.33
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.230 0.153 Inf 1.507 0.2874
## High - Low 0.523 0.140 Inf 3.746 0.0005
## Medium - Low 0.293 0.157 Inf 1.861 0.1503
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#Signifcant interaction 3.4 Growth rate G6
head(data_phenotyping)## Population Week Block ID_Rep Divsersity Population.1 Box Initial.N Start
## 1 High1 5-week Block3 52 High 1 1 30 7/8/21
## 2 High1 5-week Block3 53 High 1 2 30 7/8/21
## 3 High13 5-week Block4 129 High 13 1 30 7/15/21
## 4 High13 5-week Block4 130 High 13 2 30 7/15/21
## 5 High13 5-week Block4 131 High 13 3 30 7/15/21
## 6 High14 5-week Block4 132 High 14 1 30 7/15/21
## End Larvae Pupae Adults Lambda Genetic_Diversity Nb_ttx p_adults
## 1 8/12/21 3 9 97 3.2333333 High 109 0.8899083
## 2 8/12/21 1 1 8 0.2666667 High 10 0.8000000
## 3 8/19/21 2 0 84 2.8000000 High 86 0.9767442
## 4 8/19/21 1 2 104 3.4666667 High 107 0.9719626
## 5 8/19/21 2 1 84 2.8000000 High 87 0.9655172
## 6 8/19/21 1 0 83 2.7666667 High 84 0.9880952
## p_pupae p_larvae He_remain He_start He_exp He_end
## 1 0.08256881 0.027522936 0.9986008 0.9986008 0.9679591 0.9666047
## 2 0.10000000 0.100000000 0.9986008 0.9986008 0.9679591 0.9666047
## 3 0.00000000 0.023255814 0.9986008 0.9986008 0.9786327 0.9772634
## 4 0.01869159 0.009345794 0.9986008 0.9986008 0.9786327 0.9772634
## 5 0.01149425 0.022988506 0.9986008 0.9986008 0.9786327 0.9772634
## 6 0.00000000 0.011904762 0.9986008 0.9986008 0.9766045 0.9752381#Without considering extinct populations
#tapply(data_phenotyping$Lambda,data_phenotyping$Population,length)
#############################################################
################## Determine distribution ###################
#############################################################
## Poisson lognormal model
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
#summary(mlog)
# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
#summary(msqrt)
# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
## Compare AIC
AIC(mlog, msqrt,mnorm)## df AIC
## mlog 7 120.8153
## msqrt 7 156.0901
## mnorm 7 528.8983#Square root a better fits to the data
#But we compare different values: Growth rate and Growth rate+1
# ## Simulate data
x.teo.log <- unlist(simulate(mlog))
x.teo.sqrt <- unlist(simulate(msqrt))
x.teo.norm <- unlist(simulate(mnorm))
# ## QQplot to compared Negative binomial and Poisson log normal distributions
qqplot(data_phenotyping$Lambda, x.teo.log,
main="QQ-plot", xlab="Observed data", ylab="Simulated data",
las=1, xlim = c(0, 10), ylim = c(0, 10), col='blue') ## QQ-plot
points(sort(data_phenotyping$Lambda), sort(x.teo.sqrt), pch=1, col='red')
points(sort(data_phenotyping$Lambda), sort(x.teo.norm), pch=1, col='green')
abline(0,1)
# ## Add legend
legend(0, 7, legend=c("Log", "Sqrt", "Normal"),
col=c("blue", "red", "green"),
pch=1, bty="n")
# ## Normal distribution provides a good fit to the data
# QQ plot
# qqnorm(resid(mlog))
# qqline(resid(mlog))
# qqnorm(resid(msqrt))
# qqline(resid(msqrt))
# qqnorm(resid(mnorm))
# qqline(resid(mnorm)) #best fit
#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)
hist(residuals(msqrt),col="yellow",freq=F) #Seems the better fit
hist(residuals(mnorm),col="yellow",freq=F) #Seems the better fit
#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest ## 1-mle-norm
## "not rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest ## 1-mle-norm
## "not rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest## 1-mle-norm
## "not rejected"## Compare AIC
AIC(mlog, msqrt, mnorm)## df AIC
## mlog 7 120.8153
## msqrt 7 156.0901
## mnorm 7 528.8983#SLog a better fits to the data
#############################################################
################## Analysis ###################
#############################################################
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
summary(mod0)## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## Data: data_phenotyping
##
## REML criterion at convergence: 514.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3747 -0.6242 -0.0778 0.5409 3.4451
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2321 0.4817
## Residual 0.7480 0.8649
## Number of obs: 186, groups: Population, 57
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.569709 0.193717 13.265
## Genetic_DiversityMedium -0.873366 0.238587 -3.661
## Genetic_DiversityLow -0.700898 0.222333 -3.152
## BlockBlock4 0.220484 0.214487 1.028
## BlockBlock5 -0.003695 0.249586 -0.015
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -0.460
## Gntc_DvrstL -0.587 0.363
## BlockBlock4 -0.636 0.079 0.175
## BlockBlock5 -0.592 0.089 0.232 0.451mod1 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 5 532.53 548.66 -261.26 522.53
## mod0 7 520.90 543.48 -253.45 506.90 15.635 2 0.0004027 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 2.64 0.135 42.2 2.31 2.98
## Medium 1.77 0.198 51.5 1.28 2.26
## Low 1.94 0.177 58.2 1.51 2.38
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.873 0.239 48.0 3.653 0.0018
## High - Low 0.701 0.223 52.5 3.143 0.0076
## Medium - Low -0.172 0.261 53.0 -0.660 0.7873
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates############
###### Predict
############
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_phenotyping$Genetic_Diversity),
Block = "Block4")
data_predict_extinct$predict <- predict(mod0,
type = "response",
re.form = NA,
newdata = data_predict_extinct)3.5 Remaining He difference
Rmisc::summarySE(data[data$Generation_Eggs==2,],
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2 Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3 Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320mfull <- lm(He_remain ~ Genetic_Diversity, data=data[data$Generation_Eggs==2,])
mless <- lm(He_remain ~ 1, data=data[data$Generation_Eggs==2,])
anova(mfull,mless) #Keep effect## Analysis of Variance Table
##
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 81 0.00631
## 2 83 3.14572 -2 -3.1394 20162 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.9918 0.001698 81 0.9877 0.9960
## Medium 0.6743 0.001840 81 0.6698 0.6788
## Low 0.5404 0.001513 81 0.5367 0.5441
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.318 0.00250 81 126.829 <.0001
## High - Low 0.451 0.00227 81 198.470 <.0001
## Medium - Low 0.134 0.00238 81 56.200 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesRmisc::summarySE(data[data$Generation_Eggs==2,],
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9918493 0.001771315 0.0003408897 0.0007007089
## 2 Medium 23 0.6743040 0.006431259 0.0013410101 0.0027810847
## 3 Low 34 0.5404225 0.012690807 0.0021764555 0.0044280320#Calcul for end:
Rmisc::summarySE(data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",],
measurevar = c("He_end"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_end sd se ci
## 1 High 27 0.9663114 0.01931641 0.003717444 0.007641317
## 2 Medium 18 0.6378048 0.03493729 0.008234799 0.017373906
## 3 Low 26 0.5090656 0.03443416 0.006753094 0.013908257mfull <- lm(He_end ~ Genetic_Diversity, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
mless <- lm(He_end ~ 1, data=data[data$Generation_Eggs==2&data$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect## Analysis of Variance Table
##
## Model 1: He_end ~ Genetic_Diversity
## Model 2: He_end ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 68 0.06009
## 2 70 2.97522 -2 -2.9151 1649.3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.966 0.00572 68 0.952 0.980
## Medium 0.638 0.00701 68 0.621 0.655
## Low 0.509 0.00583 68 0.495 0.523
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.329 0.00905 68 36.316 <.0001
## High - Low 0.457 0.00817 68 55.978 <.0001
## Medium - Low 0.129 0.00912 68 14.124 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates3.6 Growth rate and He
#Best model with genetic diversity
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod1 <- lme4::lmer(Lambda ~ He_end + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod2 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod5 <- lme4::lmer(Lambda ~ exp(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
MuMIn::model.sel(mod0, mod1,mod2,mod3,mod4,mod5)## Model selection table
## (Int) Blc Gnt_Dvr He_end log(He_end) exp(He_end) family df
## mod1 0.8154 + 1.766 gaussian(identity) 6
## mod5 0.3445 + 0.8311 gaussian(identity) 6
## mod3 2.5470 + 1.259 gaussian(identity) 6
## mod0 2.5700 + + gaussian(identity) 7
## mod2 2.0770 + gaussian(identity) 5
## mod4 2.0770 + gaussian(identity) 5
## logLik AICc delta weight
## mod1 -257.506 527.5 0.00 0.408
## mod5 -258.097 528.7 1.18 0.226
## mod3 -258.156 528.8 1.30 0.213
## mod0 -257.449 529.5 2.05 0.147
## mod2 -263.551 537.4 9.95 0.003
## mod4 -263.551 537.4 9.95 0.003
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Population'#Makes sense: best model log He
x <- seq(1:100)
y <- seq(1001:1100)
plot(log(x),y)
#############################################################
################## Analysis ###################
#############################################################
mod3 <- lme4::lmer(Lambda ~ log(He_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_end),])
anova(mod3,mod4,test="Chisq") # difference## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_end) + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod4 5 532.53 548.66 -261.26 522.53
## mod3 6 522.65 542.01 -255.33 510.65 11.878 1 0.000568 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod3)## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ log(He_end) + Block + (1 | Population)
## Data: data_phenotyping[!is.na(data_phenotyping$He_end), ]
##
## REML criterion at convergence: 516.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3851 -0.6642 -0.0835 0.5336 3.4375
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2567 0.5067
## Residual 0.7473 0.8644
## Number of obs: 186, groups: Population, 57
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.547160 0.201530 12.639
## log(He_end) 1.259386 0.351126 3.587
## BlockBlock4 0.229875 0.219317 1.048
## BlockBlock5 0.004275 0.254183 0.017
##
## Correlation of Fixed Effects:
## (Intr) lg(H_) BlckB4
## log(He_end) 0.655
## BlockBlock4 -0.625 -0.154
## BlockBlock5 -0.581 -0.197 0.446#############################################################
################## Predict ###################
#############################################################
data_predict_lambda <- data.frame(He_end =
seq(min(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end),
max(data_phenotyping[!is.na(data_phenotyping$He_end),]$He_end),
0.01),
Block = "Block4")
data_predict_lambda$lambda_predict <- predict(mod3,
type = "response",
re.form = NA,
newdata = data_predict_lambda)
PLOT_He_expect <- PLOT_He_Final + geom_line(data = data_predict_lambda,
aes(x = He_end, y = lambda_predict),
linetype = "longdash", col = "grey40", size = 1.25)
# #
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"), PLOT_He_expect,
# base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)3.7 Adults G6
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 1818.1 1834.3 -904.1 1808.1 181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23572 -0.14532 0.02672 0.14825 0.41394
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.1739 0.4170
## Population (Intercept) 0.1085 0.3294
## Number of obs: 186, groups: ID_Rep:Population, 186; Population, 57
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2943 0.0812 52.884 < 2e-16 ***
## Genetic_DiversityMedium -0.5132 0.1438 -3.569 0.000358 ***
## Genetic_DiversityLow -0.3047 0.1295 -2.352 0.018660 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.563
## Gntc_DvrstL -0.625 0.356# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population/ID_Rep),
# family = "poisson", data = data_5week)
#dispersion_lme4::glmer(m0) #dispersion ok
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 1826.3 1836.0 -910.15 1820.3
## m0 5 1818.1 1834.3 -904.07 1808.1 12.164 2 0.002284 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 1818.1 1834.3 -904.1 1808.1 181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23572 -0.14532 0.02672 0.14825 0.41394
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.1739 0.4170
## Population (Intercept) 0.1085 0.3294
## Number of obs: 186, groups: ID_Rep:Population, 186; Population, 57
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2943 0.0812 52.884 < 2e-16 ***
## Genetic_DiversityMedium -0.5132 0.1438 -3.569 0.000358 ***
## Genetic_DiversityLow -0.3047 0.1295 -2.352 0.018660 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.563
## Gntc_DvrstL -0.625 0.356############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.29 0.0812 Inf 4.10 4.49
## Medium 3.78 0.1188 Inf 3.50 4.06
## Low 3.99 0.1011 Inf 3.75 4.23
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.513 0.144 Inf 3.569 0.0010
## High - Low 0.305 0.130 Inf 2.352 0.0489
## Medium - Low -0.208 0.156 Inf -1.340 0.3728
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") ## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 4 1823.6 1836.5 -907.81 1815.6
## m0 6 1820.1 1839.5 -904.07 1808.1 7.4783 2 0.02377 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") ## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 1844.6 1854.3 -919.30 1838.6
## m0 5 1838.5 1854.6 -914.26 1828.5 10.086 2 0.006455 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.8 Adults G2
m0 <- glm(Nb_adults ~ Genetic_Diversity + Block, family = "poisson", data = data[data$Generation_Eggs==2,])
summary(m0)##
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson",
## data = data[data$Generation_Eggs == 2, ])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -18.7057 -2.7621 0.2126 3.3983 11.4454
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.73890 0.01497 383.312 < 2e-16 ***
## Genetic_DiversityMedium -0.19408 0.01628 -11.920 < 2e-16 ***
## Genetic_DiversityLow -0.19278 0.01460 -13.205 < 2e-16 ***
## BlockBlock4 0.10767 0.01579 6.817 9.3e-12 ***
## BlockBlock5 0.17743 0.01600 11.088 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2379.6 on 83 degrees of freedom
## Residual deviance: 2027.9 on 79 degrees of freedom
## AIC: 2667.2
##
## Number of Fisher Scoring iterations: 4m1<- glm(Nb_adults ~ Block, family = "poisson", data = data[data$Generation_Eggs==2,])
anova(m0,m1,test="Chisq") # difference## Analysis of Deviance Table
##
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 79 2027.9
## 2 81 2241.4 -2 -213.52 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 5.834 0.01051 Inf 5.809 5.859
## Medium 5.640 0.01243 Inf 5.610 5.670
## Low 5.641 0.01022 Inf 5.617 5.666
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.1941 0.0163 Inf 11.920 <.0001
## High - Low 0.1928 0.0146 Inf 13.205 <.0001
## Medium - Low -0.0013 0.0161 Inf -0.081 0.9964
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates